function FPU_cal()
    N = 6;
    K = 1;
    omega0 = K/1;
    alpha = 0;
    t0 = 1e2; %单位为us
    step = 1e3; %演化步长
    list_len = 1e3; %总演化数
    t_min = step * t0 * 1e-6;
    t_max = step * list_len * t0 * 1e-6;
    time_list = linspace(t_min, t_max ,list_len);
    
    %计算alpha=0时的简正模式和简正频率
    [evector, omega] = eigen_cal(N);
    evector_inv = inv(evector);
    
    %初始化
    ini_state = cell(1,N);
    choice = 1; %选择：针对原坐标的初始化(0) 或 针对简正坐标的初始化(1)
    if (choice == 0)
        %针对原体系的初始化
        for i = 1:N
            ini_state{i} = [0,0];
        end
        ini_state{1} = [1,1];
    else
        %针对简正模式的初始化
        ini_qlist = zeros(1,N);
        ini_qtlist = zeros(1,N);
        ini_qlist(1) = 1;
        ini_qtlist(1) = 1;
        ini_xlist = ini_qlist * evector_inv;
        ini_xtlist = ini_qtlist * evector_inv;
        for i=1:N
            ini_state{i} = [ini_xlist(i), ini_xtlist(i)]; 
        end
    end
    en0 = tot_energy(N, ini_state, K, alpha); %计算体系的初始能量
    
    %设置输出参数
    error_list = zeros(1,list_len);
    energycell = cell(1,N); %记录各模式的能量随时间的变化
    xcell = cell(1,N); %记录各质点位移随时间的变化
    xtcell = cell(1,N); %记录各质点速度随时间的变化
    qcell = cell(1,N); %记录各简正模式位移随时间的变化
    qtcell = cell(1,N); %记录各简正模式速度随时间的变化
    for i = 1:N
        energycell{i} = zeros(1,list_len);
        xcell{i} = zeros(1,list_len);
        xtcell{i} = zeros(1,list_len);
        qcell{i} = zeros(1,list_len);
        qtcell{i} = zeros(1,list_len);
    end
    
    state = ini_state;
    
    for j = 1:list_len
        
        statej = v_verlet(N, state, t0, step, K, alpha);
        state = statej;
        
        en = tot_energy(N, statej, K, alpha);
        error_list(j) = abs(en-en0)/en0; 
        
        for k = 1:N
            xcell{k}(j) = statej{k}(1);
            xtcell{k}(j) = statej{k}(2);
        end
        
        xlist = zeros(1,N);
        xtlist = zeros(1,N);
        for k = 1:N
            xlist(k) = statej{k}(1);
            xtlist(k) = statej{k}(2);
        end
        
        for k=1:N
            qcell{k}(j) = xlist * evector(:,k);
            qtcell{k}(j) = xtlist * evector(:,k);
        end
        
        for k=1:N
            energyj = 0.5 * (qtcell{k}(j))^2 + 0.5 * ((omega(k,k)*omega0)^2) * (qcell{k}(j))^2;
            energycell{k}(j) = energyj;
        end
        
    end
    
    mode_en_error_list = zeros(1,list_len);
    for q = 1:list_len
        en_mode = 0;
        for r = 1:N
            en_mode = en_mode + energycell{r}(q);
        end
        mode_en_error_list(q) = abs(en_mode - en0)/en0;
        %mode_en_error_list(q) = en_mode;
    end
    
    %记录时间以区分不同模拟结果
    date = datestr(now, 'yyyy-mm-dd_HH-MM-SS');
    %将数据存为txt文件
    %第1个文件：时间序列、各时刻用原坐标计算出的能量与en0的相对误差、各时刻质点位移、各时刻质点速度
    fname1 = sprintf('N=%d_K=%.1f_alpha=%.3f_t0=%dus_step=%d_st=%d_ed=%d_%s-1.txt',N, K, alpha, t0, step, t_min, t_max ,date);
    fid1 = fopen(['.\',fname1],'w');
    for m = 1:list_len
        fprintf(fid1, '%d\t',time_list(m));
    end
    fprintf(fid1,'\n');
    for m = 1:list_len
        fprintf(fid1, '%d\t',error_list(m));
    end
    fprintf(fid1,'\n');
    for m=1:N
        for n=1:list_len
            fprintf(fid1, '%d\t',xcell{m}(n));
        end
        fprintf(fid1, '\n');
    end
    for m=1:N
        for n=1:list_len
            fprintf(fid1, '%d\t',xtcell{m}(n));
        end
        fprintf(fid1, '\n');
    end
    fclose(fid1);
    % % % % % % %
    %第2个文件：时间序列、各时刻用简正坐标计算出的能量（不含非线性项能量）与en0的相对误差、各时刻各模式的位移、各时刻各模式的速度
    fname2 = sprintf('N=%d_K=%.1f_alpha=%.3f_t0=%dus_step=%d_st=%d_ed=%d_%s-2.txt',N, K, alpha, t0, step, t_min, t_max, date);
    fid2 = fopen(['.\',fname2],'w');
    for m = 1:list_len
        fprintf(fid2, '%d\t',time_list(m));
    end
    fprintf(fid2,'\n');
    for m = 1:list_len
        fprintf(fid1, '%d\t',mode_en_error_list(m));
    end
    fprintf(fid1,'\n');
    for m=1:N
        for n=1:list_len
            fprintf(fid2, '%d\t',qcell{m}(n));
        end
        fprintf(fid2, '\n');
    end
    for m=1:N
        for n=1:list_len
            fprintf(fid2, '%d\t',qtcell{m}(n));
        end
        fprintf(fid2, '\n');
    end
    for m=1:N
        for n=1:list_len
            fprintf(fid2, '%d\t',energycell{m}(n));
        end
        fprintf(fid2, '\n');
    end
    fclose(fid2);
    % % % % % % %
    %第3个文件：记录结束时原坐标下的结束时间、位移和速度
    fname3 = sprintf('N=%d_K=%.1f_alpha=%.3f_t0=%dus_step=%d_endtime=%d_%s.txt',N, K, alpha, t0, step, t_max, date);
    fid3 = fopen(['.\',fname3],'w');
    for m = 1:N
        fprintf(fid3, '%d\t', t_max);
    end
    fprintf(fid3,'\n');
    for m = 1:N
        fprintf(fid3, '%d\t', xcell{m}(list_len));
    end
    fprintf(fid3,'\n');
    for m = 1:N
        fprintf(fid3, '%d\t', xtcell{m}(list_len));
    end  
    fprintf(fid3,'\n');
    fclose(fid3);    
    
    %画图
    figure(1)
    scatter(time_list, error_list,'DisplayName','x');
    hold on 
    scatter(time_list, mode_en_error_list,'DisplayName','q');
    xlabel('time(s)');
    ylabel('relative error of energy');
    legend('show');
    
%     figure(2)
%     for p=1:N
%         dotsname = sprintf('particle%d',p);
%         plot(time_list, xcell{p},'DisplayName',dotsname);
%         hold on
%     end
%     xlabel('time(s)');
%     ylabel('x(a.u.)');
%     legend('show');
%     
%     figure(3)
%     for p=1:N
%         dotsname = sprintf('particle%d',p);
%         plot(time_list, xtcell{p},'DisplayName',dotsname);
%         hold on
%     end
%     xlabel('time(s)');
%     ylabel('v(a.u.)');
%     legend('show');
    
    figure(4)
    for p=1:N
        dotsname = sprintf('mode%d',p);
        plot(time_list, qcell{p},'DisplayName',dotsname);
        hold on
    end
    xlabel('time(s)');
    ylabel('q(a.u.)');
    legend('show');
    
    figure(5)
    for p=1:N
        dotsname = sprintf('mode%d',p);
        plot(time_list, qtcell{p},'DisplayName',dotsname);
        hold on
    end
    xlabel('time(s)');
    ylabel('p(a.u.)');
    legend('show');
    
    figure(6)
    for p=1:5
        dotsname = sprintf('mode%d',p);
        plot(time_list, energycell{p},'DisplayName',dotsname);
        hold on
    end
    
    for p=96:N
        dotsname = sprintf('mode%d',p);
        plot(time_list, energycell{p},'DisplayName',dotsname);
        hold on
    end
    xlabel('time(s)');
    ylabel('energy(a.u.)');
    legend('show');

end